Robust ecological pattern formation induced by demographic noise 
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We demonstrate that demographic noise can induce persistent spatial pattern formation and 
temporal oscillations in the Levin-Segel predator-prey model for plankton-herbivore population dy- 
namics. Although the model exhibits a Turing instability in mean field theory, demographic noise 
greatly enlarges the region of parameter space where pattern formation occurs. To distinguish be- 
tween patterns generated by fluctuations and those present at the mean fleld level in real ecosystems, 
we calculate the power spectrum in the noise-driven case and predict the presence of fat tails not 
present in the mean fleld case. These results may account for the prevalence of large-scale ecological 
patterns, beyond that expected from traditional non-stochastic approaches. 

PACS numbers: 87.23.Cc, 87.10.Mn, 02.50.Ey, 05.40.-a 
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Many years ago, Turing showed how diffusion, nor- 
mally thought of as a homogenizing influence, can give 
rise to pattern- forming instabilities 1]. Only recently, 
however, have field observations provided strong support 
for the presence of Turing patterns in ecosystems, where 
diffusional processes abound, at least in principle. The 
slow moving tussock moth population in California to- 
gether with its faster moving parasites 3 as well as sev- 
eral plant-resource systems [3[ have been identified as 
satisfying, qualitatively at least, the key requirements for 
diffusion driven pattern formation. Observed patterns of 
plankton populations have also been proposed to arise 
from Turing instabilities, at least over short length scales 

mm- 

The common feature of these systems is positive feed- 
back coupled to slow diffusion (usually associated with a 
species labeled an "activator" that activates both itself 
and another species called the "inhibitor"), and nega- 
tive feedback coupled to faster diffusion associated with 
the inhibitor. This combination of diffusion and feedback 
promotes the formation of patterns, because local patches 
are promoted through positive feedback, but are only 
able to spread a limited distance before the fast diffusion 
and associated negative feedback of the inhibitor prevents 
further spread. It is hypothesized that this mechanism 
is responsible for a great deal of ecosystem level pattern 
formation 0,0]. 

One particular class of ecological pattern forming sys- 
tems, predator-prey (or organism-natural enemy) sys- 
tems has been extensively analyzed theoretically (see for 
example, t4,, li, i, i, [3 ) and is beginning to allow qual- 
itative comparison to field data along with more system 
specific theory 0, [ll| . A difficulty in directly comparing 
the results of this large body of theory to field obser- 
vations is that in many cases, models only exhibit Tur- 
ing instabilities if the predator diffusivity is much larger 
than the prey diffusivity or the parameters are fine tuned 
a i i lllj . The qualitative argument made above for 
pattern formation does not depend on very large differ- 



ences in diffusivities, nor on additional ecological details, 
and indeed, there are ecological pattern-forming systems 
which do not apparently display very large separation of 
diffusivities So what is the origin of pattern for- 

mation in such systems? 

One approach to such questions of ordering is to in- 
clude levels of detail that in some sense force the re- 
sponse of the system. For example, whereas simple mean 
field predator-prey models do not show population oscil- 
lations, they can be made to do so by the inclusion of 
predator satiation effects 12]. However, such levels of re- 
alism do not need to be invoked, because there is a sim- 
pler explanation: intrinsic or demographic noise. This 
may seem counterintuitive, because adding noise to a sys- 
tem is usually thought of as reducing ordering by adding 
entropy; and indeed, this is exactly what is observed in 
several models, such as percolation models of epidemics 
and spin models of forest canopy gaps [l^] . Surpris- 
ingly, however, systematic treatments of individual- level 
models (ILMs) of predator-prey dynamics show that the 
population fluctuations become amplified 15], and lead 
to time-dependent oscillations (quasi-cycles) that can be 
distinguished from deterministic limit cycle behavior [l6|. 
Disappointingly, to date, no novel spatial effects of de- 
mographic noise have been identified, despite several 
attempts [i3,|ll. 

In this Rapid Communication, we demonstrate that 
noise-induced pattern formation arises in a simple but 
biologically-relevant predator-prey model, and show that 
if it is analyzed as an ILM, patterns occur over a much 
larger range of ecologically relevant parameters than pre- 
dicted by MFT, even in the thermodynamic limit. We 
accomplish this by calculating the phase diagram and 
power spectrum of the model analytically. We also pre- 
dict that experimental noise driven patterns will have 
power spectra with fat tails not present in patterns driven 
by instabilities present in MFT. Finally, we show that 
quasi-cycles are also present, and we provide an interpre- 
tation of the spatiotemporal dynamics that result. 
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HEURISTIC ANALYSIS OF THE LEVIN-SEGEL 
MODEL 

Among the simplest models of ecological pattern for- 
mation was originally introduced to model plankton- 
herbivore dynamics This model takes the form 



dt^i — //V^V + blip + eV'^ — {pi + P2)i"^ 



(1) 



where the plankton population density is given by ip, 
the herbivore population density is given by (/?, hi is 
birthrate for the plankton, pi and p2 are predation, d 
is competition-driven death of the predators and e cor- 
responds to a community effect, that is the prey facili- 
tates its own birth rate. In the original presentation of 
this model, this term was intended to be a proxy for re- 
duced predator efficiency at higher prey concentrations 
It can also be interpreted as an AUee effect, wherein 
many species have enhanced reproduction at higher con- 
centrations (for a review, see [19|). From here on, we set 
Pi = Q and P2 — P for transparency of analysis. This 
does not change the qualitative results. The parame- 
ters e and d identify the prey as the activator and 
the predator as the inhibitor in the mechanism for pat- 
tern formation above and distinguish this model from 
the standard Lotka-Volterra based individual level mod- 
els recently analyzed and demonstrated not to contain 



patterns in jlTI. Il8| 



The model contains a stable homogeneous coexistence 
state when 



p > e and p'^ > de 
with fixed point populations given by 



bid 



hp 



p-^ — de' p^ 
It contains a Turing instability if 



de 



(2) 



(3) 




p/d - ^Jp/d-e/pj 



(4) 



When the model violates the stability conditions in 
Eq. [21 the plankton population diverges and a plank- 
ton regulation term (i.e. —fip^) is required to make the 
model valid. Such a term would only materially affect 
the outcomes of this analysis near the instability, where 
it would decrease the set of parameters for which pat- 
tern formation occurs. To examine the behavior of the 
model, we take the generic set of 0(1) kinetic parame- 
ters bi = 1/2, e = 1/2, d = 1/2 and p = 1. With these 



generic parameters Eq. [4] shows that non-generic diffu- 
sivities, v/fj, > 27.8, are required for pattern formation. 
Similar results are obtained for other stable, generic pa- 
rameter sets. 

Demographic noise may change this picture [2^ by in- 
hibiting the decay of transient patterns. Turing insta- 
bilities occur when, for some specific set of wave vectors, 
small perturbations no longer decay. However, we expect 
that even when the parameters are tuned away from the 
Turing instability, perturbations with some wavelengths 
may decay more slowly than others, leading to tran- 
sient patterns. Demographic noise would maintain these 
patterns by generating continual perturbations. This is 
reminiscent of extrinsic noise driven patterns reported in 
other contexts [ill, [53, [H . 



To quantify this heuristic argument, we look at the 
Fourier transformed dynamics of the fluctuations from 
the coexistence fixed point with added white noise ^, vari- 
ance 1. These dynamics are given by 

— iwx = Ax + ^ (5) 
The matrix A is the Fourier transformed stability matrix 

j^^r~iyP~p^ pip \ 

\^ —pijj — ^/c^ + eip J ^ ' 

Simple manipulations yield the power spectrum 



{pbiip + ^vk'^ — up' 



^k^ev f 1 - + u?(ie ~ p)^ -{p^ v)k^) 

\ ev J 



(7) 



Very approximately, we expect from Eq. [7] that pat- 
terns (indicated by peaks in the power spectrum) form 
whenever eu > ppL. This is much less stringent than 
Eq. [3] and can be satisfied for generic sets of parame- 
ters. However, to reliably demonstrate our hypotheses 
and extract experimental predictions, we next perform a 
systematic study of demographic noise from an individ- 
ual level model. 



INDIVIDUAL LEVEL MODEL 

We define the individual level version of the model by 
considering a locally well mixed patch of volume V . We 
consider the following reactions 



p % pp 
PH HH 



elV 

PP U ppp 
HH e H 



(8) 
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where P denotes plankton and H denotes herbivores, 
with the parameters as described above. Stochastic tra- 
jectories of H and P, enumerated by m and n respec- 
tively, are described by the master equation 



dtP(rn, n) = bi{~nP{m, n) + {n — l)P{m, n — 1)) 
+ ^[{n - l)(n - 2)P(m,n- 1) - n{n - l)P(m,n)] 

+ —{-~mnP{m, n) -f (to — \){n + l)P(m — 1, n + 1)) 

-|- — [(m -|- l)mP{m + l,n) — m{m — 1)P(to, n)] (9) 

To analyze the master equation, we map it to a path 
integral formulation of bosonic field theory and generalize 
to space [iJ, [2^ IS 27, 2§|. To add space we consider a 
lattice of patches, and random hopping for both species 
at different rates between nearest neighbor patches. The 
resulting Lagrangian density is given by 



C =xdtz + pdtp — vzSI z — ppy p — vzi^'z 
- ppiypf + 6ip(l - eP) -f - eP) 



;^2 



^zp{\-e^-P) 



v 

d o , 



(10) 



Where z, z are noise and number variables respectively 
for herbivores, and similarly, p, p are noise and number 
variables for plankton. To analyze this Lagrangian di- 
rectly is difficult, due to exponential terms and diffusive 
noise. To make progress, we derive a systematic expan- 
sion and mean field theory (MFT) in powers moti- 
vated by the fi-expansion 3i 29| . 



We assume the forms 



z 

z = Vip + Vv rj 



■'-TV 

p = Vil; + 



(11) 
(12) 



for the fields and drop terms with negative powers of y/V. 
This yields the following form of the Lagrangian 



c = Vvci + €2 + 0(1/ W) 



(13) 



Minimizing Ci in the infinite V limit yields the MFT in 
Eqs. [1] Since we've already analyzed it, we now turn to 
£2- We represent it in matrix form as 



y^dtx - y^Ax - ^Y^Y 



(14) 



The matrix A is the stability matrix we used in the 
heuristic analysis above, Eq. [Sj The matrix B is given 
by 



2pLpijj + vipk"^ 



-pifip 



(15) 



where we have Fourier transformed the equations. We 
also now note that C2 is in the form of a Lagrangian in 
the Martin-Siggia-Rose (MSR) resp onse function formal- 
ism for Langevin equations [SOl . |3l| . Thus we can extract 
coupled Langevin equations for the fiuctuations from the 
Lagrangian by applying the MSR formalism. The result- 
ing Langevin equations with the appropriate noise and 
correlations are 



— iux = Ax + 7(0;) 
Simple manipulations yield the power spectrum 



(16) 



(xixl) 



\D22?B^i-2Di2Re{D22)B2i 



ID12PP22 



\det{D)\' 



(17) 

This expression results in a rational polynomial with 
complicated coefficients that is sixth order in k in the 
numerator, and eighth in the denominator. The denom- 
inator is the same as the denominator for the heuristic 
power spectrum in Eq. [71 Alternatively, these results 
could have been obtained by a standard f2 expansion [29j 
of the master equation [51 



DISCUSSION 

Pattern formation occurs when there is a peak in 
P{k,uj) at non-zero k. This occurs if dP/dk^ > at 
fc = 0, because for large fc, the power spectrum is a de- 
creasing function and has a negative derivative. The peak 
occurs at the point where the derivative changes sign. 
Carrying out the derivative at fc = yields 



V 

- > 



p^ibp"^ + 7de) 



p e(4p4 + 5p2^g _,_ 3g;2g2 



(18) 



Eqs. [T8l[4land the stability conditions define the phase 
diagram of the model (fig. [1]). For the purposes of the 
phase diagram, we fix the parameters as above, leaving p 
and V / pas control parameters. The phase diagram shows 
that the beyond mean field corrections expand the range 
of ecologically interesting parameters in which pattern 
formation occurs greatly. 

For larger values of fc, since the denominator in Eq. [T71 
goes as the eighth power, and the numerator as the sixth 
power of fc, it is clear that 



P cx fc" 



(19) 
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This provides an experimental prediction: in regions 
II and III of the phase diagram, the power spectrum will 
have a fat tail that decays as approximately k~^. In 
region I, the power spectrum will be dominated by the 
spatially structured mean field populations, and should 
fall off much more quickly. This is analogous to the sta- 
tistical test to distinguish quasi-cycles from limit cycles 
in predator-prey populations that recently showed popu- 
lation oscillations in wolverines to be driven by finite size 
fluctuations 1^, T^ . 




We also note that the appropriate thermodynamic 
limit of the theory is not F — > oo, but rather that the 
number of patches of size V goes to infinity. Since V is 
the volume of a locally well mixed population, it should 
never be infinity for a system in which diffusion effects 
are significant. Thus the results we have presented do not 
depend on the size of the population being studied, and 
even apply to infinite populations, provided local pop- 
ulations are finite. In ecological terms, this means that 
systems in which fluctuation effects might be expected to 
be insignificant due to large populations (e.g. plankton) 
are equally likely to contain fiuctuation driven patterns 
and cycles as systems with small populations, at least 
over length scales where diffusion is a reasonable approx- 
imation for the spatial dynamics. 

The results we have given here were calculated within 
a specific model, but we expect that they will be sub- 
stantially unchanged in any model with a slow diffusing 
activator species and a faster diffusing inhibitor species. 

This work was partially supported by National Science 
Foundation grant NSF-EF-0526747. 



FIG. 1: Phase diagram over stable parameter region in p. The 
region I phase is MFT level pattern formation, the region II 
phase is noise driven pattern formation and quasi-cycles and 
region III is a spatially homogeneous phase with quasi-cycles. 

An additional feature of the model is that oscillations 
and spatial pattern formation are essentially decoupled. 
This means that the model predicts global population 
oscillations and spatial pattern formation, but not trav- 
eling waves. The mathematical origin of this can be seen 
in Eq. [T] The fc^ term with a negative coefficient at 
w = is quickly overwhelmed by the positive fc^ depen- 
dence of the Lu'^ term as the frequency begins to grow. In 
the power spectrum (fig. [J) this can be seen as the deep 
valley between the peaks in k and cu. This interpreta- 
tion is supported by preliminary simulations of an agent 
based model. 
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FIG. 2: Power spectrum with p=l, i^//i=15 
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